This is code used to make plots based on output from the “IndividualPondModels_Forecasts85” for boththe CR and YF Read in the output files for each pond/RCP combination then work on delta figures and figures like Hamlet et al. 2020
rm(list = ls())
library(tidyverse)
library(dplyr)
library(rjags)
library(ggplot2)
library(matrixStats)
set.seed(1)
Starting with the Copper River Delta
RCP 45
Mod_CRD45 <- read.csv("Corr_Mod_CRD45.csv", header = T)
# List of unique pond names from the dataframe
pond_names <- unique(Mod_CRD45$pond_name)
# Initialize an empty dataframe to store results
avgMed_2020 <- data.frame(Pond = pond_names, AvgMed_2020 = numeric(length(pond_names)))
# Loop through each pond name
for (i in seq_along(pond_names)) {
pond_name <- pond_names[i]
# Filter data for the current pond and for dates before January 1, 2021
filtered_data <- Mod_CRD45 %>%
filter(pond_name == !!pond_name, as.Date(time) < as.Date("2021-01-01"))
# Calculate the mean of 'med' (which corresponds to X50%)
avg_med <- mean(filtered_data$med, na.rm = TRUE)
# Store the result in the dataframe
avgMed_2020[i, "AvgMed_2020"] <- avg_med
}
# Print the dataframe
print(avgMed_2020)
library(lubridate)
library(dplyr)
library(tidyr)
library(ggplot2)
# List of unique pond names from the dataframe
pond_names <- unique(Mod_CRD45$pond_name)
# Initialize a dataframe to store results with years as columns
avgMed_2030 <- data.frame(Pond = pond_names)
# Loop through each year from 2030 to 2039
for (yr in 2030:2039) {
# Calculate average medians for each year and store as a new column
avgMed_col <- sapply(pond_names, function(pond_name) {
yearly_data <- Mod_CRD45 %>%
filter(pond_name == !!pond_name, year(as.Date(time)) == yr)
avg_med <- mean(yearly_data$med, na.rm = TRUE)
return(avg_med)
})
# Add the results as a new column to the dataframe
avgMed_2030[[as.character(yr)]] <- avgMed_col
}
# Print the final dataframe
print(avgMed_2030)
# Convert to long format for ggplot2
avgMed_2030_long <- avgMed_2030 %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "AvgMed")
# Plot the data
ggplot(avgMed_2030_long, aes(x = Year, y = AvgMed, color = Pond, group = Pond)) +
geom_line() +
geom_point() +
labs(title = "Average Median per Pond from 2030 to 2039",
x = "Year",
y = "Average Median") +
theme_minimal()
library(lubridate)
library(dplyr)
library(tidyr)
library(ggplot2)
# List of unique pond names from the dataframe
pond_names <- unique(Mod_CRD45$pond_name)
# Initialize a dataframe to store results with years as columns
avgMed_2060 <- data.frame(Pond = pond_names)
# Loop through each year from 2060 to 2069
for (yr in 2060:2069) {
# Calculate average medians for each year and store as a new column
avgMed_col <- sapply(pond_names, function(pond_name) {
yearly_data <- Mod_CRD45 %>%
filter(pond_name == !!pond_name, year(as.Date(time)) == yr)
avg_med <- mean(yearly_data$med, na.rm = TRUE)
return(avg_med)
})
# Add the results as a new column to the dataframe
avgMed_2060[[as.character(yr)]] <- avgMed_col
}
# Print the final dataframe
print(avgMed_2060)
# Convert to long format for ggplot2
avgMed_2060_long <- avgMed_2060 %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "AvgMed")
# Plot the data
ggplot(avgMed_2060_long, aes(x = Year, y = AvgMed, color = Pond, group = Pond)) +
geom_line() +
geom_point() +
labs(title = "Average Median per Pond from 2060 to 2069",
x = "Year",
y = "Average Median") +
theme_minimal()
library(lubridate)
library(dplyr)
library(tidyr)
library(ggplot2)
# List of unique pond names from the dataframe
pond_names <- unique(Mod_CRD45$pond_name)
# Initialize a dataframe to store results with years as columns
avgMed_2090 <- data.frame(Pond = pond_names)
# Loop through each year from 2090 to 2099
for (yr in 2090:2099) {
# Calculate average medians for each year and store as a new column
avgMed_col <- sapply(pond_names, function(pond_name) {
yearly_data <- Mod_CRD45 %>%
filter(pond_name == !!pond_name, year(as.Date(time)) == yr)
avg_med <- mean(yearly_data$med, na.rm = TRUE)
return(avg_med)
})
# Add the results as a new column to the dataframe
avgMed_2090[[as.character(yr)]] <- avgMed_col
}
# Print the final dataframe
print(avgMed_2090)
# Convert to long format for ggplot2
avgMed_2090_long <- avgMed_2090 %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "AvgMed")
# Plot the data
ggplot(avgMed_2090_long, aes(x = Year, y = AvgMed, color = Pond, group = Pond)) +
geom_line() +
geom_point() +
labs(title = "Average Median per Pond from 2090 to 2099",
x = "Year",
y = "Average Median") +
theme_minimal()
# Rename columns to reflect only the years
avgMed_2020 <- avgMed_2020 %>%
rename(AvgMed_2020 = AvgMed_2020)
avgMed_2030 <- avgMed_2030 %>%
rename(`2030` = `2030`)
avgMed_2060 <- avgMed_2060 %>%
rename(`2060` = `2060`)
avgMed_2090 <- avgMed_2090 %>%
rename(`2090` = `2090`)
# Combine the data frames
combined_avgMed <- avgMed_2020 %>%
left_join(avgMed_2030, by = "Pond") %>%
left_join(avgMed_2060, by = "Pond") %>%
left_join(avgMed_2090, by = "Pond")
# Print the final dataframe
print(combined_avgMed)
# Convert to long format for ggplot2, excluding AvgMed_2020
avgMed_long <- combined_avgMed %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "AvgMed") %>%
filter(Year %in% as.character(2021:2199))
# Plot the data
ggplot(avgMed_long, aes(x = Year, y = AvgMed, color = Pond, group = Pond)) +
geom_line() +
geom_point() +
labs(title = "Average Median per Pond in the CRD - RCP 4.5",
x = "Year",
y = "Average Median") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),# Rotate x-axis labels by 45 degrees
plot.title = element_text(hjust = 0.5))
library(dplyr)
# Assuming combined_avgMed is loaded and structure verified
Delta_combined_avgMed <- combined_avgMed # Use your actual dataframe
# Columns to calculate differences for
years <- c("2030", "2031", "2032", "2033", "2034",
"2035", "2036", "2037", "2038", "2039",
"2060", "2061", "2062", "2063", "2064",
"2065", "2066", "2067", "2068", "2069",
"2090", "2091", "2092", "2093", "2094",
"2095", "2096", "2097", "2098", "2099")
# Loop through each year and calculate the difference
for (year in years) {
# Create a new column for the difference
diff_col_name <- paste0("diff_", year)
Delta_combined_avgMed[[diff_col_name]] <- Delta_combined_avgMed[[year]] - Delta_combined_avgMed$`AvgMed_2020`
}
# Display the structure of Delta_combined_avgMed to verify
str(Delta_combined_avgMed)
'data.frame': 9 obs. of 62 variables:
$ Pond : chr "BVS" "CAB" "TIN" "SQR" ...
$ AvgMed_2020: num 5.65 5.89 7.64 8.4 8.45 ...
$ 2030 : num 7 7.24 9 9.76 9.84 ...
$ 2031 : num 7.98 8.21 9.97 10.75 10.81 ...
$ 2032 : num 9.19 9.44 11.2 11.97 12.03 ...
$ 2033 : num 7.49 7.73 9.49 10.27 10.32 ...
$ 2034 : num 8.21 8.44 10.2 10.97 11.03 ...
$ 2035 : num 4.74 4.98 6.73 7.52 7.58 ...
$ 2036 : num 6.1 6.34 8.08 8.9 8.93 ...
$ 2037 : num 8.6 8.83 10.6 11.35 11.42 ...
$ 2038 : num 7.41 7.67 9.41 10.18 10.24 ...
$ 2039 : num 8.16 8.41 10.16 10.92 10.97 ...
$ 2060 : num 9.49 9.73 11.49 12.27 12.32 ...
$ 2061 : num 9.33 9.57 11.32 12.08 12.16 ...
$ 2062 : num 9.9 10.1 11.9 12.7 12.7 ...
$ 2063 : num 10.8 11.1 12.8 13.6 13.6 ...
$ 2064 : num 8.14 8.4 10.16 10.92 10.97 ...
$ 2065 : num 10.2 10.4 12.2 13 13 ...
$ 2066 : num 9.03 9.27 11.03 11.8 11.86 ...
$ 2067 : num 8.53 8.77 10.52 11.3 11.36 ...
$ 2068 : num 9.23 9.46 11.21 11.98 12.05 ...
$ 2069 : num 10.1 10.4 12.1 12.9 12.9 ...
$ 2090 : num 9.64 9.9 11.66 12.42 12.48 ...
$ 2091 : num 11.2 11.5 13.2 14 14 ...
$ 2092 : num 11.8 12.1 13.8 14.6 14.7 ...
$ 2093 : num 10.6 10.8 12.6 13.4 13.4 ...
$ 2094 : num 10.9 11.2 12.9 13.7 13.7 ...
$ 2095 : num 10.9 11.1 12.9 13.7 13.7 ...
$ 2096 : num 11.5 11.7 13.5 14.3 14.3 ...
$ 2097 : num 10.6 10.8 12.6 13.4 13.4 ...
$ 2098 : num 10.5 10.8 12.5 13.3 13.4 ...
$ 2099 : num 11.6 11.8 13.6 14.4 14.4 ...
$ diff_2030 : num 1.35 1.35 1.37 1.36 1.38 ...
$ diff_2031 : num 2.32 2.32 2.33 2.34 2.36 ...
$ diff_2032 : num 3.54 3.54 3.57 3.57 3.57 ...
$ diff_2033 : num 1.84 1.84 1.86 1.86 1.87 ...
$ diff_2034 : num 2.56 2.54 2.56 2.56 2.57 ...
$ diff_2035 : num -0.914 -0.919 -0.903 -0.888 -0.878 ...
$ diff_2036 : num 0.444 0.449 0.446 0.495 0.477 ...
$ diff_2037 : num 2.94 2.94 2.97 2.94 2.97 ...
$ diff_2038 : num 1.75 1.77 1.78 1.77 1.78 ...
$ diff_2039 : num 2.51 2.52 2.52 2.52 2.52 ...
$ diff_2060 : num 3.84 3.84 3.85 3.86 3.87 ...
$ diff_2061 : num 3.67 3.67 3.69 3.68 3.7 ...
$ diff_2062 : num 4.24 4.24 4.26 4.27 4.27 ...
$ diff_2063 : num 5.16 5.18 5.18 5.17 5.18 ...
$ diff_2064 : num 2.49 2.5 2.52 2.52 2.51 ...
$ diff_2065 : num 4.53 4.54 4.54 4.55 4.56 ...
$ diff_2066 : num 3.38 3.37 3.4 3.39 3.4 ...
$ diff_2067 : num 2.88 2.87 2.89 2.9 2.91 ...
$ diff_2068 : num 3.57 3.57 3.57 3.58 3.59 ...
$ diff_2069 : num 4.45 4.46 4.47 4.48 4.48 ...
$ diff_2090 : num 3.99 4 4.02 4.01 4.03 ...
$ diff_2091 : num 5.55 5.57 5.56 5.56 5.59 ...
$ diff_2092 : num 6.18 6.18 6.19 6.2 6.22 ...
$ diff_2093 : num 4.93 4.93 4.95 4.96 4.95 ...
$ diff_2094 : num 5.27 5.29 5.29 5.3 5.29 ...
$ diff_2095 : num 5.24 5.24 5.26 5.27 5.26 ...
$ diff_2096 : num 5.84 5.84 5.86 5.85 5.87 ...
$ diff_2097 : num 4.93 4.94 4.95 4.96 4.96 ...
$ diff_2098 : num 4.87 4.86 4.89 4.89 4.9 ...
$ diff_2099 : num 5.94 5.95 5.96 5.96 5.97 ...
# Create a new dataframe with the calculated means and standard deviations
Delta_CRD45 <- Delta_combined_avgMed %>%
mutate(
Mean_2030_2039 = rowMeans(select(., starts_with("diff_203")), na.rm = TRUE), # Calculate row-wise mean for columns starting with "diff_203"
SD_2030_2039 = rowSds(as.matrix(select(., starts_with("diff_203"))), na.rm = TRUE), # Calculate row-wise SD for columns starting with "diff_203"
Mean_2060_2069 = rowMeans(select(., starts_with("diff_206")), na.rm = TRUE), # Calculate row-wise mean for columns starting with "diff_206"
SD_2060_2069 = rowSds(as.matrix(select(., starts_with("diff_206"))), na.rm = TRUE), # Calculate row-wise SD for columns starting with "diff_206"
Mean_2090_2099 = rowMeans(select(., starts_with("diff_209")), na.rm = TRUE), # Calculate row-wise mean for columns starting with "diff_209"
SD_2090_2099 = rowSds(as.matrix(select(., starts_with("diff_209"))), na.rm = TRUE) # Calculate row-wise SD for columns starting with "diff_209"
)
# Print the structure of the new dataframe
str(Delta_CRD45)
'data.frame': 9 obs. of 68 variables:
$ Pond : chr "BVS" "CAB" "TIN" "SQR" ...
$ AvgMed_2020 : num 5.65 5.89 7.64 8.4 8.45 ...
$ 2030 : num 7 7.24 9 9.76 9.84 ...
$ 2031 : num 7.98 8.21 9.97 10.75 10.81 ...
$ 2032 : num 9.19 9.44 11.2 11.97 12.03 ...
$ 2033 : num 7.49 7.73 9.49 10.27 10.32 ...
$ 2034 : num 8.21 8.44 10.2 10.97 11.03 ...
$ 2035 : num 4.74 4.98 6.73 7.52 7.58 ...
$ 2036 : num 6.1 6.34 8.08 8.9 8.93 ...
$ 2037 : num 8.6 8.83 10.6 11.35 11.42 ...
$ 2038 : num 7.41 7.67 9.41 10.18 10.24 ...
$ 2039 : num 8.16 8.41 10.16 10.92 10.97 ...
$ 2060 : num 9.49 9.73 11.49 12.27 12.32 ...
$ 2061 : num 9.33 9.57 11.32 12.08 12.16 ...
$ 2062 : num 9.9 10.1 11.9 12.7 12.7 ...
$ 2063 : num 10.8 11.1 12.8 13.6 13.6 ...
$ 2064 : num 8.14 8.4 10.16 10.92 10.97 ...
$ 2065 : num 10.2 10.4 12.2 13 13 ...
$ 2066 : num 9.03 9.27 11.03 11.8 11.86 ...
$ 2067 : num 8.53 8.77 10.52 11.3 11.36 ...
$ 2068 : num 9.23 9.46 11.21 11.98 12.05 ...
$ 2069 : num 10.1 10.4 12.1 12.9 12.9 ...
$ 2090 : num 9.64 9.9 11.66 12.42 12.48 ...
$ 2091 : num 11.2 11.5 13.2 14 14 ...
$ 2092 : num 11.8 12.1 13.8 14.6 14.7 ...
$ 2093 : num 10.6 10.8 12.6 13.4 13.4 ...
$ 2094 : num 10.9 11.2 12.9 13.7 13.7 ...
$ 2095 : num 10.9 11.1 12.9 13.7 13.7 ...
$ 2096 : num 11.5 11.7 13.5 14.3 14.3 ...
$ 2097 : num 10.6 10.8 12.6 13.4 13.4 ...
$ 2098 : num 10.5 10.8 12.5 13.3 13.4 ...
$ 2099 : num 11.6 11.8 13.6 14.4 14.4 ...
$ diff_2030 : num 1.35 1.35 1.37 1.36 1.38 ...
$ diff_2031 : num 2.32 2.32 2.33 2.34 2.36 ...
$ diff_2032 : num 3.54 3.54 3.57 3.57 3.57 ...
$ diff_2033 : num 1.84 1.84 1.86 1.86 1.87 ...
$ diff_2034 : num 2.56 2.54 2.56 2.56 2.57 ...
$ diff_2035 : num -0.914 -0.919 -0.903 -0.888 -0.878 ...
$ diff_2036 : num 0.444 0.449 0.446 0.495 0.477 ...
$ diff_2037 : num 2.94 2.94 2.97 2.94 2.97 ...
$ diff_2038 : num 1.75 1.77 1.78 1.77 1.78 ...
$ diff_2039 : num 2.51 2.52 2.52 2.52 2.52 ...
$ diff_2060 : num 3.84 3.84 3.85 3.86 3.87 ...
$ diff_2061 : num 3.67 3.67 3.69 3.68 3.7 ...
$ diff_2062 : num 4.24 4.24 4.26 4.27 4.27 ...
$ diff_2063 : num 5.16 5.18 5.18 5.17 5.18 ...
$ diff_2064 : num 2.49 2.5 2.52 2.52 2.51 ...
$ diff_2065 : num 4.53 4.54 4.54 4.55 4.56 ...
$ diff_2066 : num 3.38 3.37 3.4 3.39 3.4 ...
$ diff_2067 : num 2.88 2.87 2.89 2.9 2.91 ...
$ diff_2068 : num 3.57 3.57 3.57 3.58 3.59 ...
$ diff_2069 : num 4.45 4.46 4.47 4.48 4.48 ...
$ diff_2090 : num 3.99 4 4.02 4.01 4.03 ...
$ diff_2091 : num 5.55 5.57 5.56 5.56 5.59 ...
$ diff_2092 : num 6.18 6.18 6.19 6.2 6.22 ...
$ diff_2093 : num 4.93 4.93 4.95 4.96 4.95 ...
$ diff_2094 : num 5.27 5.29 5.29 5.3 5.29 ...
$ diff_2095 : num 5.24 5.24 5.26 5.27 5.26 ...
$ diff_2096 : num 5.84 5.84 5.86 5.85 5.87 ...
$ diff_2097 : num 4.93 4.94 4.95 4.96 4.96 ...
$ diff_2098 : num 4.87 4.86 4.89 4.89 4.9 ...
$ diff_2099 : num 5.94 5.95 5.96 5.96 5.97 ...
$ Mean_2030_2039: num 1.83 1.83 1.85 1.85 1.86 ...
$ SD_2030_2039 : num 1.29 1.29 1.3 1.29 1.29 ...
$ Mean_2060_2069: num 3.82 3.82 3.84 3.84 3.85 ...
$ SD_2060_2069 : num 0.803 0.807 0.802 0.803 0.805 ...
$ Mean_2090_2099: num 5.27 5.28 5.29 5.3 5.3 ...
$ SD_2090_2099 : num 0.642 0.641 0.637 0.639 0.643 ...
# Optionally, view the first few rows to check the new columns
head(Delta_CRD45)
# Load necessary libraries
library(dplyr)
# Select and rename columns from Delta_CRD45
CRD45_decadal <- Delta_CRD45 %>%
select(Pond, starts_with("Mean_"), starts_with("SD_"))
# Rename columns without the 'Y' prefix for means
names(CRD45_decadal) <- sub("^Mean_", "Mean_", names(CRD45_decadal))
# Rename columns without the 'Y' prefix for standard deviations
names(CRD45_decadal) <- sub("^SD_", "SD_", names(CRD45_decadal))
# Print the structure of the updated dataframe
str(CRD45_decadal)
'data.frame': 9 obs. of 7 variables:
$ Pond : chr "BVS" "CAB" "TIN" "SQR" ...
$ Mean_2030_2039: num 1.83 1.83 1.85 1.85 1.86 ...
$ Mean_2060_2069: num 3.82 3.82 3.84 3.84 3.85 ...
$ Mean_2090_2099: num 5.27 5.28 5.29 5.3 5.3 ...
$ SD_2030_2039 : num 1.29 1.29 1.3 1.29 1.29 ...
$ SD_2060_2069 : num 0.803 0.807 0.802 0.803 0.805 ...
$ SD_2090_2099 : num 0.642 0.641 0.637 0.639 0.643 ...
# Optionally, view the first few rows to check the new columns
head(CRD45_decadal)
# Load required packages
library(ggplot2)
library(gridExtra)
library(ggpubr)
# Create a grouped bar plot with error bars for Mean 2030
plot_2030 <- ggplot(CRD45_decadal, aes(x = Pond, y = Mean_2030_2039)) +
geom_bar(stat = "identity", fill = "#FFD92F", color = "black") +
geom_errorbar(aes(ymin = Mean_2030_2039 - SD_2030_2039, ymax = Mean_2030_2039 + SD_2030_2039),
width = 0.4, # Adjust the width of error bars as needed
position = position_dodge(width = 0.9)) +
labs(x = "", y = "Delta Temp (C)") +
ylim(-1, 5) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(plot.title = element_text(hjust = 0.5)) # Centered title
plot_2030
# Create a grouped bar plot with error bars for Mean 2060
plot_2060 <- ggplot(CRD45_decadal, aes(x = Pond, y = Mean_2060_2069)) +
geom_bar(stat = "identity", fill = "#F46D43", color = "black") +
geom_errorbar(aes(ymin = Mean_2060_2069 - SD_2060_2069, ymax = Mean_2060_2069 + SD_2060_2069),
width = 0.4, # Adjust the width of error bars as needed
position = position_dodge(width = 0.9)) +
labs(x = "", y = "") + # No y-axis label
ylim(-1, 5) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(plot.title = element_text(hjust = 0.5)) + # Centered title
theme(axis.text.y = element_blank(), # Remove Y-axis tick mark labels
axis.ticks.y = element_blank()) # Remove Y-axis ticks
plot_2060
# Create a grouped bar plot with error bars for Mean 2090
plot_2090 <- ggplot(CRD45_decadal, aes(x = Pond, y = Mean_2090_2099)) +
geom_bar(stat = "identity", fill = "#D73027", color = "black") +
geom_errorbar(aes(ymin = Mean_2090_2099 - SD_2090_2099, ymax = Mean_2090_2099 + SD_2090_2099),
width = 0.4, # Adjust the width of error bars as needed
position = position_dodge(width = 0.9)) +
labs(x = "", y = "") + # No y-axis label
ylim(-1, 5) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(plot.title = element_text(hjust = 0.5)) + # Centered title
theme(axis.text.y = element_blank(), # Remove Y-axis tick mark labels
axis.ticks.y = element_blank()) # Remove Y-axis ticks
plot_2090
Warning: Removed 9 rows containing missing values or values outside the scale range (`geom_bar()`).
# Combine plots into a panel plot
Delta_CRD45_plot <- ggarrange(plot_2030, plot_2060, plot_2090, ncol = 3)
Warning: Removed 9 rows containing missing values or values outside the scale range (`geom_bar()`).
# Save the combined plot
ggsave("Corr_DeltaCRD45.png", plot = Delta_CRD45_plot, width = 8, height = 3)
Delta_CRD45
# Load required libraries
library(dplyr)
library(broom)
library(tidyr)
# Step 1: Prepare the data
data_long <- Delta_CRD45 %>%
select(Pond, matches("^diff_203[0-9]$")) %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "Difference")
boxplot(data_long$Difference ~ data_long$Pond)
# Step 2: Conduct ANOVA by grouping by Pond
res_CRD45 <- aov(Difference ~ Pond, data = data_long)
summary(res_CRD45)
Df Sum Sq Mean Sq F value Pr(>F)
Pond 8 0.01 0.0012 0.001 1
Residuals 81 135.53 1.6732
# Step 3: Tukey HSD
post_test_CRD45 <- TukeyHSD(res_CRD45)
post_test_CRD45
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Difference ~ Pond, data = data_long)
$Pond
diff lwr upr p adj
CAB-BVS 1.214803e-03 -1.842514 1.844943 1
CAN-BVS 2.805931e-02 -1.815669 1.871788 1
EYS-BVS 2.289556e-02 -1.820833 1.866624 1
RHM-BVS 1.206862e-03 -1.842522 1.844935 1
SQR-BVS 2.018862e-02 -1.823540 1.863917 1
TIN-BVS 1.588325e-02 -1.827845 1.859612 1
TIS-BVS 1.889326e-02 -1.824835 1.862622 1
WDD-BVS 2.287431e-02 -1.820854 1.866603 1
CAN-CAB 2.684451e-02 -1.816884 1.870573 1
EYS-CAB 2.168076e-02 -1.822048 1.865409 1
RHM-CAB -7.940660e-06 -1.843737 1.843721 1
SQR-CAB 1.897381e-02 -1.824755 1.862702 1
TIN-CAB 1.466845e-02 -1.829060 1.858397 1
TIS-CAB 1.767845e-02 -1.826050 1.861407 1
WDD-CAB 2.165951e-02 -1.822069 1.865388 1
EYS-CAN -5.163752e-03 -1.848892 1.838565 1
RHM-CAN -2.685245e-02 -1.870581 1.816876 1
SQR-CAN -7.870694e-03 -1.851599 1.835858 1
TIN-CAN -1.217606e-02 -1.855905 1.831553 1
TIS-CAN -9.166053e-03 -1.852895 1.834563 1
WDD-CAN -5.185002e-03 -1.848914 1.838544 1
RHM-EYS -2.168870e-02 -1.865417 1.822040 1
SQR-EYS -2.706942e-03 -1.846436 1.841022 1
TIN-EYS -7.012306e-03 -1.850741 1.836716 1
TIS-EYS -4.002301e-03 -1.847731 1.839726 1
WDD-EYS -2.124991e-05 -1.843750 1.843707 1
SQR-RHM 1.898175e-02 -1.824747 1.862710 1
TIN-RHM 1.467639e-02 -1.829052 1.858405 1
TIS-RHM 1.768640e-02 -1.826042 1.861415 1
WDD-RHM 2.166745e-02 -1.822061 1.865396 1
TIN-SQR -4.305364e-03 -1.848034 1.839423 1
TIS-SQR -1.295359e-03 -1.845024 1.842433 1
WDD-SQR 2.685692e-03 -1.841043 1.846414 1
TIS-TIN 3.010005e-03 -1.840719 1.846739 1
WDD-TIN 6.991056e-03 -1.836738 1.850720 1
WDD-TIS 3.981051e-03 -1.839748 1.847710 1
Delta_CRD45
# Load required libraries
library(dplyr)
library(broom)
# Step 1: Prepare the data
data_long <- Delta_CRD45 %>%
select(Pond, matches("^diff_206[0-9]$")) %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "Difference")
boxplot(data_long$Difference ~ data_long$Pond)
# Step 2: Conduct ANOVA by grouping by Pond
res_CRD45 <- aov(Difference ~ Pond, data = data_long)
summary(res_CRD45)
Df Sum Sq Mean Sq F value Pr(>F)
Pond 8 0.01 0.0013 0.002 1
Residuals 81 52.61 0.6495
# Step 3: Tukey HSD
post_test_CRD45 <- TukeyHSD(res_CRD45)
post_test_CRD45
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Difference ~ Pond, data = data_long)
$Pond
diff lwr upr p adj
CAB-BVS 0.0035901603 -1.145125 1.152306 1
CAN-BVS 0.0283876339 -1.120328 1.177103 1
EYS-BVS 0.0195832180 -1.129132 1.168299 1
RHM-BVS -0.0013213523 -1.150037 1.147394 1
SQR-BVS 0.0188165617 -1.129899 1.167532 1
TIN-BVS 0.0160150233 -1.132701 1.164731 1
TIS-BVS 0.0177892348 -1.130926 1.166505 1
WDD-BVS 0.0292115139 -1.119504 1.177927 1
CAN-CAB 0.0247974736 -1.123918 1.173513 1
EYS-CAB 0.0159930576 -1.132723 1.164709 1
RHM-CAB -0.0049115127 -1.153627 1.143804 1
SQR-CAB 0.0152264014 -1.133489 1.163942 1
TIN-CAB 0.0124248630 -1.136291 1.161140 1
TIS-CAB 0.0141990745 -1.134517 1.162915 1
WDD-CAB 0.0256213535 -1.123094 1.174337 1
EYS-CAN -0.0088044159 -1.157520 1.139911 1
RHM-CAN -0.0297089862 -1.178425 1.119007 1
SQR-CAN -0.0095710722 -1.158287 1.139145 1
TIN-CAN -0.0123726106 -1.161088 1.136343 1
TIS-CAN -0.0105983991 -1.159314 1.138117 1
WDD-CAN 0.0008238800 -1.147892 1.149539 1
RHM-EYS -0.0209045703 -1.169620 1.127811 1
SQR-EYS -0.0007666563 -1.149482 1.147949 1
TIN-EYS -0.0035681947 -1.152284 1.145147 1
TIS-EYS -0.0017939831 -1.150510 1.146922 1
WDD-EYS 0.0096282959 -1.139087 1.158344 1
SQR-RHM 0.0201379140 -1.128578 1.168854 1
TIN-RHM 0.0173363756 -1.131379 1.166052 1
TIS-RHM 0.0191105872 -1.129605 1.167826 1
WDD-RHM 0.0305328662 -1.118183 1.179248 1
TIN-SQR -0.0028015384 -1.151517 1.145914 1
TIS-SQR -0.0010273269 -1.149743 1.147688 1
WDD-SQR 0.0103949522 -1.138321 1.159111 1
TIS-TIN 0.0017742115 -1.146941 1.150490 1
WDD-TIN 0.0131964906 -1.135519 1.161912 1
WDD-TIS 0.0114222790 -1.137293 1.160138 1
Delta_CRD45
# Load required libraries
library(dplyr)
library(broom)
# Step 1: Prepare the data
data_long <- Delta_CRD45 %>%
select(Pond, matches("^diff_209[0-9]$")) %>%
pivot_longer(cols = -Pond, names_to = "Year", values_to = "Difference")
boxplot(data_long$Difference ~ data_long$Pond)
# Step 2: Conduct ANOVA by grouping by Pond
res_CRD45 <- aov(Difference ~ Pond, data = data_long)
summary(res_CRD45)
Df Sum Sq Mean Sq F value Pr(>F)
Pond 8 0.01 0.0014 0.003 1
Residuals 81 33.39 0.4122
# Step 3: Tukey HSD
post_test_CRD45 <- TukeyHSD(res_CRD45)
post_test_CRD45
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Difference ~ Pond, data = data_long)
$Pond
diff lwr upr p adj
CAB-BVS 0.005808524 -0.9092851 0.9209022 1
CAN-BVS 0.031775890 -0.8833178 0.9468695 1
EYS-BVS 0.020810076 -0.8942836 0.9359037 1
RHM-BVS -0.001850173 -0.9169438 0.9132435 1
SQR-BVS 0.023209929 -0.8918837 0.9383036 1
TIN-BVS 0.018623321 -0.8964703 0.9337170 1
TIS-BVS 0.012662910 -0.9024307 0.9277566 1
WDD-BVS 0.024980879 -0.8901128 0.9400745 1
CAN-CAB 0.025967366 -0.8891263 0.9410610 1
EYS-CAB 0.015001551 -0.9000921 0.9300952 1
RHM-CAB -0.007658697 -0.9227524 0.9074350 1
SQR-CAB 0.017401405 -0.8976922 0.9324951 1
TIN-CAB 0.012814797 -0.9022789 0.9279085 1
TIS-CAB 0.006854386 -0.9082393 0.9219480 1
WDD-CAB 0.019172354 -0.8959213 0.9342660 1
EYS-CAN -0.010965815 -0.9260595 0.9041278 1
RHM-CAN -0.033626063 -0.9487197 0.8814676 1
SQR-CAN -0.008565961 -0.9236596 0.9065277 1
TIN-CAN -0.013152570 -0.9282462 0.9019411 1
TIS-CAN -0.019112980 -0.9342066 0.8959807 1
WDD-CAN -0.006795012 -0.9218887 0.9082986 1
RHM-EYS -0.022660248 -0.9377539 0.8924334 1
SQR-EYS 0.002399854 -0.9126938 0.9174935 1
TIN-EYS -0.002186755 -0.9172804 0.9129069 1
TIS-EYS -0.008147165 -0.9232408 0.9069465 1
WDD-EYS 0.004170803 -0.9109229 0.9192645 1
SQR-RHM 0.025060102 -0.8900336 0.9401538 1
TIN-RHM 0.020473494 -0.8946202 0.9355671 1
TIS-RHM 0.014513083 -0.9005806 0.9296067 1
WDD-RHM 0.026831051 -0.8882626 0.9419247 1
TIN-SQR -0.004586609 -0.9196803 0.9105070 1
TIS-SQR -0.010547019 -0.9256407 0.9045466 1
WDD-SQR 0.001770949 -0.9133227 0.9168646 1
TIS-TIN -0.005960411 -0.9210541 0.9091332 1
WDD-TIN 0.006357558 -0.9087361 0.9214512 1
WDD-TIS 0.012317968 -0.9027757 0.9274116 1
# Extract unique pond names from the dataframe
short_names <- unique(Mod_CRD45$pond_name)
# Initialize an empty dataframe to store monthly mean results
monthly_means <- data.frame(Pond = character(), Month = character(), AvgTemp = numeric(), Year = numeric(), stringsAsFactors = FALSE)
# Convert 'time' column to Date class
Mod_CRD45$time <- as.Date(Mod_CRD45$time, format="%Y-%m-%d")
# Extract year and month from 'time' column
Mod_CRD45$Year <- format(Mod_CRD45$time, "%Y")
Mod_CRD45$Month <- format(Mod_CRD45$time, "%B")
# Loop through each unique pond name
for (short_name in short_names) {
# Filter data for the current pond
pond_data <- subset(Mod_CRD45, pond_name == short_name)
# Loop through each month
for (month in month.name) {
# Filter data for the current month
filtered_data <- subset(pond_data, Month == month)
# Calculate the mean of 'med' for the current month
avg_med <- mean(filtered_data$med, na.rm = TRUE)
# Append results to monthly_means dataframe
monthly_means <- rbind(monthly_means, data.frame(Pond = short_name, Month = month, AvgTemp = avg_med, Year = 2020))
}
}
# Print the monthly means dataframe
print(monthly_means)
# Plotting the data
ggplot(monthly_means, aes(x = Month, y = AvgTemp, color = Pond, group = Pond)) +
geom_line() +
geom_point() +
scale_x_discrete(limits = month.name) +
labs(title = "Monthly Average Temperatures for Each Pond",
x = "Month",
y = "Average Temperature",
color = "Pond") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# Define the correct order for months
month_order <- c("January", "February", "March", "April", "May", "June",
"July", "August", "September", "October", "November", "December")
# Recalculate monthly averages with correct month ordering
calculate_monthly_averages <- function(data, decade_start, decade_end) {
# Initialize empty dataframe to store monthly averages
monthly_avg_table <- data.frame(Pond = character(), Month = factor(), AvgTemp = numeric(), Year = integer(), stringsAsFactors = FALSE)
# Loop through each year in the decade
for (yr in decade_start:decade_end) {
# Filter data for the current year
yearly_data <- data %>%
filter(year(time) == yr)
# Loop through each pond
for (pond_name in unique(data$pond_name)) {
pond_data <- yearly_data %>%
filter(pond_name == pond_name)
# Loop through each month
for (month in month_order) {
# Filter data for the current month
monthly_data <- pond_data %>%
filter(format(time, "%B") == month)
# Calculate the mean of 'med' for the current month
avg_temp <- mean(monthly_data$med, na.rm = TRUE)
# Append results to monthly_avg_table dataframe
monthly_avg_table <- rbind(monthly_avg_table, data.frame(Pond = pond_name, Month = factor(month, levels = month_order), AvgTemp = avg_temp, Year = yr))
}
}
}
return(monthly_avg_table)
}
# Calculate monthly averages for each decade
monthly_avgX50_table_2030s <- calculate_monthly_averages(Mod_CRD45, 2030, 2039)
monthly_avgX50_table_2060s <- calculate_monthly_averages(Mod_CRD45, 2060, 2069)
monthly_avgX50_table_2090s <- calculate_monthly_averages(Mod_CRD45, 2090, 2099)
# Add a Decade column to each dataframe
monthly_means <- monthly_means %>%
mutate(Decade = "2020s")
monthly_avgX50_table_2030s <- monthly_avgX50_table_2030s %>%
mutate(Decade = "2030s")
monthly_avgX50_table_2060s <- monthly_avgX50_table_2060s %>%
mutate(Decade = "2060s")
monthly_avgX50_table_2090s <- monthly_avgX50_table_2090s %>%
mutate(Decade = "2090s")
# Combine all monthly average tables into one dataframe
combined_monthly_avg <- bind_rows(
monthly_means,
monthly_avgX50_table_2030s,
monthly_avgX50_table_2060s,
monthly_avgX50_table_2090s
)
# Print the combined dataframe
print(combined_monthly_avg)
# Define the reference year
reference_year <- 2020
# Check if the column exists and is correctly populated
summary(combined_monthly_avg)
Pond Month AvgTemp Year Decade
Length:3348 Length:3348 Min. :-1.329 Min. :2020 Length:3348
Class :character Class :character 1st Qu.: 7.762 1st Qu.:2036 Class :character
Mode :character Mode :character Median :11.285 Median :2064 Mode :character
Mean :10.913 Mean :2063
3rd Qu.:14.112 3rd Qu.:2092
Max. :17.872 Max. :2099
# Calculate average temperature for each month in the reference year
ref_year_avg <- combined_monthly_avg %>%
filter(Year == reference_year) %>%
group_by(Month) %>%
summarise(avg_temp_ref_year = mean(AvgTemp, na.rm = TRUE))
# Join reference year averages with the main dataset
combined_monthly_avg_with_diff <- combined_monthly_avg %>%
left_join(ref_year_avg, by = "Month") %>%
mutate(diff_AvgTemp = AvgTemp - avg_temp_ref_year) %>%
select(-avg_temp_ref_year) # Remove the reference year average column
# Display the structure to verify
str(combined_monthly_avg_with_diff)
'data.frame': 3348 obs. of 6 variables:
$ Pond : chr "BVS" "BVS" "BVS" "BVS" ...
$ Month : chr "January" "February" "March" "April" ...
$ AvgTemp : num 4.55 3.92 3.45 5.65 8.25 ...
$ Year : num 2020 2020 2020 2020 2020 2020 2020 2020 2020 2020 ...
$ Decade : chr "2020s" "2020s" "2020s" "2020s" ...
$ diff_AvgTemp: num -1.58 -1.59 -1.6 -1.65 -1.74 ...
# Filter out data for the reference year
DeltaMonthly_CRD45 <- combined_monthly_avg_with_diff %>%
filter(Year != reference_year)
# Display the filtered data
print(DeltaMonthly_CRD45)
NA
library(dplyr)
DeltaMonthly_CRD45
DiffMonthly_CRD45 <- DeltaMonthly_CRD45 %>%
select(Pond, Month, Year, diff_AvgTemp)
# Create a new column based on year ranges
DiffMonthly_CRD45$Decade <- ifelse(DiffMonthly_CRD45$Year >= 2030 & DiffMonthly_CRD45$Year < 2040, "2030s",
ifelse(DiffMonthly_CRD45$Year >= 2060 & DiffMonthly_CRD45$Year < 2070, "2060s",
ifelse(DiffMonthly_CRD45$Year >= 2090 & DiffMonthly_CRD45$Year < 2100, "2090s",
"Other")))
# Print the updated data frame
print(DiffMonthly_CRD45)
# Convert Month to Date type
DiffMonthly_CRD45 <- DiffMonthly_CRD45 %>%
mutate(Month = as.Date(paste0(Month, " 1 ", Year), format = "%B %d %Y"))
print(DiffMonthly_CRD45)
# Group by Pond and Month, then calculate mean and standard deviation of diff_AvgTemp
summary_stats <- DiffMonthly_CRD45 %>%
dplyr::mutate(Month = format(Month, "%m")) %>% # Extracts YYYY-MM format
dplyr::group_by(Pond, Month, Decade) %>%
dplyr::summarise(
mean_diff_AvgTemp = mean(diff_AvgTemp, na.rm = TRUE),
sd_diff_AvgTemp = sd(diff_AvgTemp, na.rm = TRUE)
)
`summarise()` has grouped output by 'Pond', 'Month'. You can override using the `.groups` argument.
# View the summary statistics
print(summary_stats)
# Group by just Month, then calculate mean and standard deviation of diff_AvgTemp
summary_stats2 <- DiffMonthly_CRD45 %>%
dplyr::mutate(Month = format(Month, "%m")) %>% # Extracts YYYY-MM format
dplyr::group_by(Month, Decade) %>%
dplyr::summarise(
mean_diff_AvgTemp = mean(diff_AvgTemp, na.rm = TRUE),
sd_diff_AvgTemp = sd(diff_AvgTemp, na.rm = TRUE)
)
`summarise()` has grouped output by 'Month'. You can override using the `.groups` argument.
# View the summary statistics
print(summary_stats2)
Plotting the monthly temperature change graphs
2030s
summary_stats # grouped by month and pond
summary_stats2 # grouped by just month
# 2. Monthly averages across ponds
MonthlyStats_2030 <- summary_stats2 %>%
filter(Decade == "2030s")
ggplot(data = MonthlyStats_2030, aes(x = Month)) +
geom_point(aes(y = mean_diff_AvgTemp)) +
theme_classic()
# Convert Month from character to numeric to ensure proper ordering in ggplot
MonthlyStats_2030$Month <- as.numeric(MonthlyStats_2030$Month)
# Plot using ggplot
MonthlyStats_2030Plot <- ggplot(data = MonthlyStats_2030, aes(x = Month)) +
geom_line(aes(y = mean_diff_AvgTemp)) + # Point plot with mean values
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Horizontal dashed line at y = 0
geom_ribbon(aes(ymin = mean_diff_AvgTemp - sd_diff_AvgTemp, ymax = mean_diff_AvgTemp + sd_diff_AvgTemp), fill = "blue", alpha = 0.3) + # Ribbon for positive and negative error based on sd values
scale_x_continuous(breaks = 1:12, labels = month.abb) + # Adjust x-axis labels to show month abbreviations
ylim(-4,6)+
theme_classic() +
labs(x = "", y = "Delta Temp (C)") +
theme(axis.text.x = element_text(angle = 65, hjust = 1)) +
theme(axis.title = element_text(size = 16),
axis.text = element_text(size = 14))
MonthlyStats_2030Plot
2060s
# 2. Monthly averages across ponds
MonthlyStats_2060 <- summary_stats2 %>%
filter(Decade == "2060s")
ggplot(data = MonthlyStats_2060, aes(x = Month)) +
geom_point(aes(y = mean_diff_AvgTemp)) +
theme_classic()
# Convert Month from character to numeric to ensure proper ordering in ggplot
MonthlyStats_2060$Month <- as.numeric(MonthlyStats_2060$Month)
# Plot using ggplot
MonthlyStats_2060Plot <- ggplot(data = MonthlyStats_2060, aes(x = Month)) +
geom_line(aes(y = mean_diff_AvgTemp)) + # Point plot with mean values
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Horizontal dashed line at y = 0
geom_ribbon(aes(ymin = mean_diff_AvgTemp - sd_diff_AvgTemp, ymax = mean_diff_AvgTemp + sd_diff_AvgTemp), fill = "blue", alpha = 0.3) + # Ribbon for positive and negative error based on sd values
scale_x_continuous(breaks = 1:12, labels = month.abb) + # Adjust x-axis labels to show month abbreviations
ylim(-4,6)+
theme_classic() +
labs(x = "", y = "") +
theme(axis.text.x = element_text(angle = 65, hjust = 1)) +
theme(axis.title = element_text(size = 16),
axis.text = element_text(size = 14))+
theme(axis.text.y = element_blank(), # Remove Y-axis tick mark labels
axis.ticks.y = element_blank())
MonthlyStats_2060Plot
2090s
# 2. Monthly averages across ponds
MonthlyStats_2090 <- summary_stats2 %>%
filter(Decade == "2090s")
ggplot(data = MonthlyStats_2090, aes(x = Month)) +
geom_point(aes(y = mean_diff_AvgTemp)) +
theme_classic()
# Convert Month from character to numeric to ensure proper ordering in ggplot
MonthlyStats_2090$Month <- as.numeric(MonthlyStats_2090$Month)
# Plot using ggplot
MonthlyStats_2090Plot <- ggplot(data = MonthlyStats_2090, aes(x = Month)) +
geom_line(aes(y = mean_diff_AvgTemp)) + # Point plot with mean values
geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Horizontal dashed line at y = 0
geom_ribbon(aes(ymin = mean_diff_AvgTemp - sd_diff_AvgTemp, ymax = mean_diff_AvgTemp + sd_diff_AvgTemp), fill = "blue", alpha = 0.3) + # Ribbon for positive and negative error based on sd values
scale_x_continuous(breaks = 1:12, labels = month.abb) + # Adjust x-axis labels to show month abbreviations
ylim(-4,6)+
theme_classic() +
theme(axis.text.x = element_text(angle = 65, hjust = 1)) +
labs(x = "", y = "") +
theme(axis.title = element_text(size = 16),
axis.text = element_text(size = 14))+
theme(axis.text.y = element_blank(), # Remove Y-axis tick mark labels
axis.ticks.y = element_blank())
MonthlyStats_2090Plot
Combine the plots together here for now
# Combine plots into a panel plot
DeltaMonth_CRD45_plot <- ggarrange(MonthlyStats_2030Plot, MonthlyStats_2060Plot, MonthlyStats_2090Plot, ncol = 3)
# Save the combined plot
ggsave("Corr_DeltaMonthCRD45.png", plot = DeltaMonth_CRD45_plot, width = 12, height = 3)
Create .csv files of the data for these graphs
write.csv(MonthlyStats_2030, file = "Corr_MonthlyStatsCRD45_2030.csv")
write.csv(MonthlyStats_2060, file = "Corr_MonthlyStatsCRD45_2060.csv")
write.csv(MonthlyStats_2090, file = "Corr_MonthlyStatsCRD45_2090.csv")